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We present a strong constraint on variation of the proton-to-electron mass ratio, /i, over cosmological 
time scales using molecular hydrogen transitions in optical quasar spectra. Using high quality spectra 
of quasars Q0405— 443, Q0347— 383 and Q0528— 250, variation in ^ relative to the present day value is 
limited to A/i//i — (2.6 ± 3.0) x 10^®. We reduce systematic errors compared to previous works by 
substantially improving the spectral wavelength calibration method and by fitting absorption profiles to 
the forest of hydrogen Lyman a transitions surrounding each H2 transition. Our results are consistent with 
no variation, and inconsistent with a previous ^ 4cr detection of /i variation involving Q0405— 443 and 
Q0347— 383. If the results of this work and those suggesting that a may be varying are both correct, then 
this would tend to disfavour certain grand unification models. 

PACS numbers: 98.62.Ra, 95.30.Dr, 98.80.Es, 14.20.Dh 
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Searches have been undertaken in recent years for cos- 
mological variations in fundamental, dimensionless con- 
stants. These searches are motivated by predictions of 
Kaluza-Klein theory, string theory and other grand unifi- 
cation theories that the so-called "fundamental constants" 
may evolve over cosmological time scales. Although much 
of the focus has been on a, the fine structure constant, 
others have examined the proton-to-electron mass ratio, 
fi = mp/rrie. The quantum chromodynamical scale, to 
which fi is sensitive, may vary faster than the quantum 
electrodynamical scale, hence ^ may vary more than a ||T1. 
The wavelengths of the Lyman and Werner transitions of 
the H2 molecule are sensitive to ^u, and examination of H2 
absorption systems in quasar spectra allows one to search 
for any such variation, as was first noted by [2J. 

Attempts from 1995 to 2004 to detect a variation in 
IJL yielded results statistically consistent with no change 
Il3lll|5l|6l|71[8l|9l[ll[nl- These searches (with the excep- 
tion of Ivanchik et al. ifTTI ) were impaired by insufficiently 
accurate laboratory measurements of the H2 wavelengths, 
as well as lower quality quasar spectra. Recent laboratory 
XUV laser measurements lfT2l[T3]| have yielded substantial 
improvements in H2 wavelength accuracy. 



Using these newly available wavelengths, Reinhold et al. 
( [2006) |14| reanalysed the observed H2 wavelengths de- 
rived by Ivanchik et al. lITTl from Very Large Telecope 
(VLT) spectra of absorbers associated with Q0405— 443 (at 
redshift z w 2.595) and Q0347-383 (at z w 3.025), find- 
ing a change in ^ of Afi/fj, = (2.4 it 0.6) x 10^^, where 
A/i//i = (^2 — /Uo)/^o> /^z is the measured value of fi 
at redshift z, and /ip is the present day laboratory value. 
However, it has since been demonstrated lITSl that the tech- 
niques used to calibrate the wavelength scale of UVES 



(Ultraviolet and Visual Echelle Spectrograph, on the VLT) 
produce both long and short range calibration errors [15|. 
These calibration errors directly impact the calculation of 
A^/fi. It is therefore important to re-analyse these spec- 
tra using the improved wavelength calibration techniques 
of 1 15 1, and we do so here. We also analyse an absorber to- 
wards Q0528— 250 (at z ~ 2.811), which provides a new, 
strong constraint on Afi/fi. We use the Voigt profile fitting 
program VPFIT to analyse our spectra. 

For a given H2 transition observed in an absorbing cloud 
at redshift Zabs: the first-order shift in the wavelength A^ 
compared to the laboratory wavelength Ag is given by 



A, = Ao(l + 2abs)(l + i^zA/x/^) 



(1) 



where Ki is the sensitivity coefficient associated with each 
transition, given by Ki = (din Xi) / {din fx). Zabs is the 
redshift of the transitions measured provided that Afi/fi = 
0. If Afi/fi 7^ 0, Zahs corresponds to the redshift, de- 
termined from the ensemble of available transitions, of 
a transitions with Ki = 0. Previous works have cal- 
culated Ki within a semi-empirical framework llT4l [T6l ; 
Reinhold et al. [14J recently produced Ki coefficients of 
improved accuracy by including effects beyond the Bom- 
Oppenheimer Approximation. We use the Ki coefficients 
calculated by Reinhold et al. f\A\ and Ubachs et al. lfT6ll . 

For a series of H2 transitions, the best-fit value of Ajjl/ jjl 
may be determined in one of two ways. In the first, the 
"reduced redshift method" ("RRM"), one calculates an ob- 
served redshift Zi for each transition, and then defines a 
reduced redshift 







^abs 



1 + ^abs 



A/i 



K,. 



(2) 



A/i//i can then be determined by a linear fit to the ob- 
served C,i vs Ki distribution. 

The second method ("VPFIT method") involves fitting 
all available transitions simultaneously and solving for a 
single redshift for each identifiable absorbing H2 compo- 
nent in the system. We refer to each fitted redshift as 
a velocity component because of the close proximity of 
these components in velocity space. Afi/fi is estimated 
by perturbing the laboratory H2 wavelengths as Aq -^ 
Ao(l + Ki^An/fi)) and minimising x^ for the spectral 
data fitted. The value of Afi/fi at the minimum x^ is the 
best-fit value. 

Although the VPFIT method has previously [61 been 
used to construct x^ vs Afi/fi curves, from which the 
best value of A^/fi can be estimated, we instead include 
Afi/fi as a free parameter in the fit (within VPFIT), to 
be solved for concomitantly with the other line param- 
eters; this yields a substantial improvement in computa- 
tional speed and robustness. 

The RRM has been used in most previous measurements 
and was the method used by [14|. This method is appeal- 
ing because the required numerical methods are relatively 
simple. However, the VPFIT method is preferable in that 
fewer parameters are required to fit the data. In particular, 
the VPFIT method has n„(nt — 1) fewer free parameters, 
where riy is the number of velocity components and rif is 
the number of transitions used. For Q0405— 443 the VP- 
FIT method yields 51 fewer parameters, for Q0347— 373 it 
yields 67 fewer parameters, and for Q0528— 250 it yields 
252 fewer parameters. 

It should be noted that the VPFIT method also improves 
the stability of the fitting process. In systems with multiple 
velocity components, particular transitions may have very 
poorly constrained line parameters, despite the fact that the 
best-fit line parameters may be well constrained over many 
transitions. In the RRM, this can cause the fitting algorithm 
to reject certain velocity components in some transitions, 
rendering those transitions unsuitable for inclusion in the 
fit. Using the VPFIT method, the reduction in the num- 
ber of free parameters as a result of requiring each velocity 
component to occur at a single redshift helps to stabilise 
the fit, allowing for the inclusion of a greater number of 
transitions. 

Each of the molecular hydrogen transitions involved falls 
within the Lyman a forest, a dense series of absorption 
lines blueward of the hydrogen Lyman a emission line of 
the quasar. These transitions substantially complicate the 
analysis; the narrow molecular hydrogen absorption lines 
are often situated deep within much broader, and usually 
complex, Lyman a lines. These contaminating atomic Ly- 
man a transitions are insensitive to a change in fi. In con- 
trast to previous works, we fit absorption profiles to all of 
the Lyman a transitions in the vicinity of each H2 transi- 
tion. This allows the inclusion of a greater number of H2 
transitions which would otherwise be excluded from the fit. 

We have reduced the total number of free parameters in 



the fit by tying together physically linked parameters. In 
particular, we tie the Doppler (linewidth) parameters to- 
gether for H2 transitions with the same rotational quantum 
number J of the initial state. For the transitions we have 
analysed, J € [0, 4]. Each error estimate is multipled by 
^/xl (where xt is the x^ per degree of freedom for the 
whole fit), to account for a non-ideal fit. 

Our inclusion of the Lyman a forest within the fitting 
process increases the number of free parameters in the fit 
substantially (to over 1000 in each quasar spectrum). With 
such a large parameter space, convergence of the optimiza- 
tion algorithm must be checked. [17, 18] demonstrated that 
the results of l|T9ll20ll2Tll22]| (in relation to a) were flawed 
in this respect. Our algorithm demonstrates the proper con- 
vergence, in that A/i/fi vs x^ curves possess the correct 
parabolic shape, with derived la error bounds on Afi/fi 
that agree with those produced by VPFIT. 

The system toward Q0347 — 373 contains a single H2 ve- 
locity component. The system towards Q0405— 443 has a 
second velocity component, separated by ~ 13 kms^^ in 
velocity space llT4l . However, many of its transitions are 
weak or are heavily blended, and so we have not utilised 
the second component here. 

The H2 system towards Q0528— 250 is more compli- 
cated. Previous attempts to examine this system have 
yielded varied and comparatively poor results [6^, TO'l be- 
cause the spectra used had substantially lower signal-to- 
noise ratio than those currently available from VLT/UVES. 
Ledoux et al. |23 1 report the detection of multiple velocity 
components in the Q0528— 250 system (that is, multiple 
systems separated in velocity space), and Srianand et al. 
[24| model the absorber with 2 components. 

We have tried modeling the absorber towards 
Q0528— 250 with 2, 3 and 4 velocity components. Two 
velocity components are plainly obvious as a substantial 
asymmetry in every line. Using the F test, the probability 
that the reduction in x^ from using three components 
instead of two is due to chance is p = 4 x 10^^^. That is, 
a three component model is very strongly preferred to a 
two component model. 

Comparing a model with four velocity components to a 
model with three gives p = 1.8 x 10^^. That is, the four 
component model is preferred over both the two and three 
component models. The use of a five component model 
produces a fit that is highly unstable numerically, and so 
we use the four component model as the fiducial model. 

In modeling the multiple velocity components, we re- 
quire that the ratio of the column densities between each 
velocity component is the same for transitions with the 
same quantum number J. Certain line parameters for 
the four-parameter fit (with appropriate redshifts tied) are 
given in table 11] 

Our results are set out in tables HH and HiH 

We re-sampled the d vs Ki graph with the bootstrap 
method [25|, to check for consistency. That is, we gen- 
erated 10^ new data sets by randomly drawing data points. 



Component Relative column density (cm 


'^) Redshift 


Quasar spectrum A^//i — RRM 


'2 

Xu z 


n 


1 0.00 


2.8110036(24) 


Q0405-443 (10.9 ± 7.1) x 10"" 


1.01 2.595 


52 


2 -0.10 ±0.03 


2.8111229(15) 


Q0347-373 (6.4 ± 10.3) x 10~" 


1.13 3.025 


68 


3 -0.60 ±0.10 


2.8109334(37) 


Q0405 + Q0347 (8.5 ± 5.7) x 10"^ 


1.06 2.810 


120 


4 -1.89 ±0.76 


2.8112139(91) 


Q0528-250 n/a 


n/a n/a 


64 



TABLE I: Relative column densities and redshifts for the best 
four— component fit to Q0528— 250, with la uncertainties, de- 
rived from the J = 1 set of lines. The relative column density 
is the difference between the logarithm of the column density for 
each component and the logarithm of the column density for the 
strongest component. We use the J = 1 set of lines because they 
are the largest fraction of the data set. 



Quasar spectrum 


A/i/^ — VPFIT 


xi z 


n 


Q0405-443 


(10.1 ±6.2) x 10"'' 


1.42 2.595 


52 


Q0347-373 


(8.2 ±7.4) X 10~" 


1.28 3.025 


68 


Q0528-250 


(-1.4 ±3.9) X 10"*' 


1.22 2.811 


64 


Weighted mean 


(2.6 ±3.0) X 10"" 


n/a 2.81 


n/a 



TABLE IL Values of A^/^ obtained using the VPFIT method, 
derived from each of the quasar spectra, n is the number of tran- 
sitions (per velocity component, for Q0528— 250). The weighted 
mean given here is our preferred result. 



with replacement, from the original data set, such that each 
of the new data sets has the same number of points as the 
original data set. We then obtained the slope of the linear 
fit to each of these data sets, giving A/i//i for each set; 
the mean of this ensemble should be consistent with the 
generating data. For each absorber, we found consistency 
between the bootstrap method, the VPFIT method and the 
RRM. 

The RRM is not appropriate for Q0528— 250 because 
the line parameters for each of the velocity components 
within a given transition are strongly correlated. So, for 
Q0528— 250, groups of 4 points in a Qi vs Ki plot (the 
RRM method) are not independent. Thus the RRM method 
breaks down, as it uses linear least squares fitting, which 
requires independence of all the data points. This demon- 
strates the superiority of the VPFIT method over the RRM. 
These correlations are correctly incorporated into the VP- 
FIT method, as A^//i is determined from x^ on the total 
model. 

We checked that using 4 velocity components (instead of 
2 or 3) has no significant effect on the result by solving for 
A/i//i (within VPFIT), with 2 and 3 velocity component 
models. This produces A/i//i = (—0.6 it 3.8) x 10~" and 
A/i//i = (-1.4 ± 3.8) X 10"" respectively. These are 
similar to each other and to the result for the four velocity 
component model. This demonstrates the insensitivity of 
the result to having used the statistically preferred velocity 
structure. 

Combining the three measurements of A/x//i obtained 
within VPFIT using a weighted mean yields the value 
A/i//i = (2.6 ± 3.0) X 10"". This is null at a la confi- 
dence level. This is our main result from a combined anal- 



TABLE III: Values of A/i//i using the RRM, derived from each 
of the quasar spectra, n is the number of transitions. 



ysis of all three quasar absorbers. 

For comparison with Reinhold et al. lfT4l . a reduced 
redshift plot (Fig. [T) including only Q0405-443 and 
Q0347-373 produces the result A^//i = (8.5 ± 5.7) x 
10"" (weighted fit) and A/i/^ = (7.9 ±8.1) x 10"" 
(unweighted fit). We also attempt to compare with lfT4l 
by including in the fit only those transitions used in that 
paper. For Q0405— 443, this removes 16 transitions and 
adds 3, the latter of which appear to be contaminated 
and which were excluded from our main analysis. This 
yields, from the RRM, a Q0405-443 result of A////U = 
(10.2 ± 8.9) X 10"". This is offset from the result of 
(27.8 ± 8.8) X 10"^ in El. For Q0347-373, we remove 
35 transitions that are not used in lfT4l . and include 4 which 
appear to be contaminated, to give a result of A;u//i = 
(12.0 ±14.0) X 10"", compared with (20.6 ±7.9) x 10"" 
from lfT4ll . The weighted mean of our results in this cir- 
cumstance is A/i//i = (10.7 ± 7.5) X 10"". It is diffi- 
cult to make a direct statistical comparison, due to the fact 
that the spectra analysed are not independent, however in 
both cases we see a shift of A/i//z towards 0. Although 
the inclusion of Q0528— 250 clearly shifts the combined 
Q0405-443 -i- Q0347-373 result towards zero, our com- 
bined Q0405-443 -i- Q0347-373 result is null under all 
the circumstances considered. 

Figure s|2] and |3] show reduced redshift plots not included 
in the shorter, published version of this paper. Example fits 
also not in the shorter version of this paper can be found in 
figsg|5l|6]|7l[8]|9]and[l0l 

Our final result of A^/^ = (2.6 ± 3.0) x 10"" rep- 
resents a significant increase in precision over previous 
works (a factor of ~ 2). This result is entirely consis- 
tent with A/i/// = over cosmological timescales. It is 
also consistent with the recently-published work of Mur- 
phy et al. EH, who find that A^/^ = (0.74 ± 0.47stat ± 
0.76sys) X 10"" using the inversion transitions of ammonia. 
Note, however, that the ammonia constraint is at z = 0.685 
while all our constraints are at z > 2.5; they may not be di- 
rectly compared without a theory of cosmologically evolv- 
ing ^. 

The unification of all interactions clearly requires that 
any cosmological variations in the various fundamental 
constants will be linked to each other. Grand Unified the- 
ories typically predict A///^ ps RAa/a ll27l |28l |29l, 
where both the sign and magnitude of R are strongly model 
dependent. \R\ ~ 30 — 40 emerges from many GUT mod- 



■K^ O 




FIG. 1: Reduced redshift plot (0 vs Ki) for Q0405-443 and Q0347-373 with gradient A^/fi = (8.5 ± 5.7) x 10"^ (dashed line). 
Q0347— 373 is represented by closed circles, and Q0405— 443 is represented by open circles. The unweighted fit, a dot-dashed line, is 
obscured by the weighted fit. The dotted lines give the la confidence limit on the regression line. Note that this is not our preferred 
result, because it does not include Q0528— 250 in the fit. This graph may be compared directly with fig. 2 of |il4i. 



els ||22l|28l|29|. Generally speaking, \R\ > 1. The most 
reliable constraint on a variation at present is Aa/a = 
(-5.7ifc 1.01 xlO~^ 11301; the works of lfT9l 1201 lin 1221 and 
others have been demonstrated to be unreliable ifTTl [T8l . 
Taking both this and our new null result at face value, any 
variation in /j, is almost 2 orders of magnitude below that 
expected on the basis of the a-variation results. If both 
these results are both correct, those Grand Unified models 
which predict |i?| ^ 1 are disfavoured. 

We would like to acknowledge V. V. Flambaum, for ad- 
vice and support in various stages of this work. 
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FIG. 2: Reduced redshift {Q vs Ki) for Q0405-443. This graph has gradient Afi/^i = (1.09 ± 0.71) x 10"^ (weighted fit). The 
weighted, Hne of best fit is given by the dashed line with the unweighted Une of best given by the dot-dashed line. The dotted hnes 
indicate the la (weighted) confidence intervals on the regression line. This figure is not included in the shorter version of this paper. 
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FIG. 3: Reduced redshift {Q vs Ki) for Q0347-373. This graph has gradient Afi/fi = (6.4± 10.3) x 10"*^ (weighted fit). The weighted, 
line of best fit is given by the dashed line with the unweighted line of best given by the dot-dashed line. The dotted lines indicate the 
1(7 (weighted) confidence intervals on the regression line. This figure is not included in the shorter version of this paper. 
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FIG. 4: Fitting region 12 of 45 for Q0347— 373. The centre positions of transitions are indicated by the blue tick marks. The top plot 
indicates the residuals, normalised to la (the horizontal lines). In this case, the second tick mark from the left (^ 4226A) is a H2 
transition. All lines within the fitting region except those indicated are fitted as hydrogen Lyman a. This figure is not included in the 
shorter version of this paper. 
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FIG. 5: Fitting region 24 of 45 for Q0347— 373. The centre positions of transitions are indicated by the blue tick marks. The top plot 
indicates the residuals, normalised to la (the horizontal lines). In this case, the sharp transitions located at ^ 4058. 6A and ^ 4059.1 A 
are H2 transitions. All lines within the fitting region except those indicated are fitted as hydrogen Lyman a. This figure is not included 
in the shorter version of this paper. 
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FIG. 6: Fitting region 14 of 45 for Q0405-'443. The centre positions of transitions are indicated by the blue tick marks. The top 
plot indicates the residuals, normalised to la (the horizontal lines). In this case, the sharp transition located at ^ 3779. 85A is a H2 
transition. All lines within the fitting region except those indicated are fitted as hydrogen Lyman a. This figure is not included in the 
shorter version of this paper. 
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FIG. 7: Fitting region 16 of 45 for Q0405— 443. The centre positions of transitions are indicated by the blue tick marks. The top plot 
indicates the residuals, normalised to la (the horizontal lines). In this case, the very prominent, sharp transitions located at ^ 3733. SA 
and « 3739. 9A are H2 transitions. All lines within the fitting region except those indicated are fitted as hydrogen Lyman a. This figure 
is not included in the shorter version of this paper. 
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FIG. 8: Fitting region 11 of 49 for Q0528— 250. The centre positions of transitions are indicated by the blue tick marks. The top plot 
indicates the residuals, normalised to la (the horizontal lines). In this case, the four tick tick marks from the left indicate the four 
H2 velocity components of Q0528— 250, with the right-most tick mark being a Lyman a transition. This figure is not included in the 
shorter version of this paper. 
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FIG. 9: Fitting region 24 of 49 for Q0528— 250. The centre positions of transitions are indicated by the blue tick marks. The top plot 
indicates the residuals, normalised to la (the horizontal lines). In this case, there are two sets of H2 transitions, with four velocity 
components each, clustered about ^ 3764. 4A and ^ 3765. 2A. There is an additional set of H2 transitions at ^ 3763. 2A, however these 
are not included within this fitting region, but rather in the adjacent one. All lines within the fitting region except those indicated are 
fitted as hydrogen Lyman a. This figure is not included in the shorter version of this paper. 
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FIG. 10: A comparison of the results of this work, with that of Reinhold et al. [141 and Murphy et al. 



